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Abstract. The glass transition of supercooled fluids is a particular challenge for com- 
puter simulation, because the (longest) relaxation times increase by about 15 decades 
upon approaching the transition temperature T g . Brute- force molecular dynamics sim- 
ulations, as presented here for molten Si02 and coarse-grained bead-spring models of 
polymer chains, can yield very useful insight about the first few decades of this slowing 
down. Hence this allows to access the temperature range around T c of the so-called 
mode coupling theory, whereas the dynamics around the experimental glass transi- 
tion is completely out of reach. While methods such as "parallel tempering" improve 
the situation somewhat, a method that allows to span a significant part of the region 
T g < T < T c is still lacking. Only for abstract models such as the infinite range 10-state 
Potts glass with a few hundred spins this region can be explored. However this model 
suffers from very strong finite size effects thus making it difficult to extrapolate the 
results obtained for the finite system sizes to the thermodynamic limit. 

For the case of polymer melts, two different strategies to use lattice models instead 
of continuum models are discussed: In the first approach, a mapping of an atomistically 
realistic model of polyethylene to the bond fluctuation model with suitable effective 
potentials and a temperature-dependent time rescaling factor is attempted. In the 
second approach, devoted to a test of the entropy theory, moves that are artificial but 
which lead to a faster relaxation ("slithering snake" algorithm) are used, to get at 
least static properties at somewhat lower temperatures than possible with a "realistic" 
dynamics. The merits and shortcomings of all these approaches are discussed. 



1 Introduction 

The reason for the slowing down of the dynamics of supercooled liquids and the 
resulting glass transition to an amorphous solid is one of the biggest unsolved 
problems in the physics of condensed matter]! and it is also a particular 

challenge for computer simulation [l^fl^ fiol 1 l| , |l 2| . The present introductory 
section intends to remind the reader on the main experimental facts and some 
theoretical ideas about the glass transition, and will also serve to make clear 
why in this problem there exists a gap of time-scales that simulations need to 
bridge. 
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As is well known, it is already a problem to characterize the static structure 
of a glass: the structure of an amorphous material is not regular like a crystalline 
solid, but shows only short range order similar to a liquid. However, the latter 
flows, while the amorphous solid is rigid! In fact, if one makes a scattering 
experiment, it is hard to distinguish from the structure whether one has a fluid 
above the glass transition temperature T g or a solid below T g (Fig. @) [l§|l§ . 
If one approaches T g , the structural relaxation time r which is related to the 
viscosity rj(T), for instance - increases smoothly by up to 15 decades, as shown 
schematically in Fig. |^, without any accompanying significant structural change 
detectable by scattering experiments (Fig. |l|). This increase of r](T) is often fitted 
to the Vogel-Fulcher relation jjj 

r](T) cx exp[E VF /(T - T VF )\ , (1) 

where Eyp is an effective activation barrier. From this functional form it is 
clear that rj(T) is predicted to diverge at the Vogel-Fulcher temperature Tvf, 
which is lower than T g , of course, if one invokes the empirical definition of T g via 
•q(T = T g ) = 10 13 Poise However, it is questionable whether the temperature 
dependence of 77 as given by Eq. (jl]) really holds. 

Another very common relation often used to describe various relaxation 
functions of glassforming fluids is the Kohlrausch- Williams- Watts function, also 
called stretched exponential, jjj, 

<p(t) cx exp[-(t/r)^ . (2) 

This relation involves an exponent /3 < 1, whose precise physical significance 
is somewhat obscure. Again it is unknown under which circumstances (if any) 
Eq. (^|) is exact, and in which it is just a convenient fitting formula to represent 
data. 

Often it is claimed that the glass transition is a purely kinetic phenomenon, 
and if one would be able to wait long enough (which could mean times like the 
age of the universe, however!) one could see that glass is not really a solid but still 
a fluid that flows. However, this idea is not generally accepted, since there are 
some indications that there may be an underlying quasi-equilibrium phase tran- 
sition between metastable phases, namely from the supercooled metastable fluid 
to an (also metastable) ideal glass phase (the stable phase for temperatures lower 
than the melting temperature T m is of course the crystal). Such an indication is 
Kauzmann's entropy paradox [fL5| : By studying the difference in entropy between 
liquid and crystal one finds that near T g the difference AS(T) = Sfl u id — ^crystal 
has decreased to about 1/3 of its value S m at the melting/crystallization temper- 
ature T m . If this trend is extrapolated (linearly in T) to even lower temperatures, 
AS(T) would become negative below the Kauzmann temperature To (which is 
usually quite close to the Vogel-Fulcher-temperature Tvf), see Fig. |[ It would 
indeed be paradox if the entropy of the supercooled fluid (with its disordered 
structure) were less than the entropy of the ordered solid! One possibility to 
bypass the problem is to assume that this "entropy catastrophe" is avoided by a 
phase transition at Tq (or at some temperature in between Tq and T g ). In fact, 
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for the glass transition of polymer melts Gibbs and Di Marzio p/fl] proposed 
an approximate theory that shows such a vanishing of the entropy at To, and 
subsequently Adam and Gibbs Q suggested arguments to show that Eq. (Q) 
holds with Typ — Tq. However, although these concepts enjoy some popularity, 
all these arguments are based on very crude and hardly justifiable assumptions 
and approximations, and hence they are not accepted by many researchers. For 
instance, the mode coupling theory of the glass transition (MCT)|^] claims that 
there is indeed an underlying transition but this is not a phase transition in the 
sense of thermodynamics but rather a "dynamical transition" from an ergodic 
to a nonergodic behavior. This transition should occur at a critical temperature 
T c and can be seen in the form of the time dependence of the correlation func- 
tion of density fluctuations or its Fourier transform <P q {T), see Fig. | Above T c , 
this correlator decays to zero, but as T c is approached a plateau develops whose 
"lifetime" gets larger and larger until it diverges, in the ideal case: the system 
gets "stuck", the decay of stops at the "nonergodicity parameter" f c , an 

order parameter for the glass transition that appears discontinuously at T c . The 
physical idea behind this theory is the "cage picture" : the motion of any atom 
in a dense fluid is constrained by its neighbors, which form a cage around it. At 
low enough temperatures the escape out of the cage gets blocked. MCT predicts 
that close to this dynamical transition r and r)(T) show a power-law divergence 
as one approaches T c , 

t oc r,{T) oc (T — T c )^ . (3) 

In reality this dependence is, however, observed only in a limited temperature 
interval. The way out of this dilemma is the argument that one must not neglect 
(as "idealized" mode coupling theory does ||) thermally activated processes, 
so-called "hopping processes" , by which atoms supposedly can escape from their 
cage when T is less than T c . The theory then claims Jl9| that a simple Arrhenius 
behavior results in this region, logr oc 1/T for T < T c , and in the vicinity of T c 
the power-law divergence of Eq. (|^) is rounded off to a smooth crossover from 
the power law to the Arrhenius divergence. Thus this theory does not involve 
any phase transition, there is just a smooth crossover from one type of dynamical 
behavior to another one near T c , and T g means that relaxation times have grown 
so large that the system falls out of equilibrium. 

In real systems this crossover seems to occur at a viscosity somewhere be- 
tween 10 and 10 3 Poise, i.e. a time window that can be explored with molecular 
dynamics (MD) simulations. Hence such simulations are able to investigate the 
beginning of the approach to the critical temperature T c which MCT describes 
H and hence are very useful to check the validity of this theory. However, the 
following 10 decades of the viscosity between T c and T g are out of reach for MD 
simulations so far. Unfortunately this is precisely the region that one needs to 
explore, for a definite distinction between the theories! 

Thus although straightforward atomistic MD methods @,|l| clearly face a 
dilemma, we shall nevertheless describe how far one can push this type of ap- 
proach, choosing Si02 as an example (Sec. 2). A method to extend the range of 
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times and accessible temperatures somewhat is the concept of "parallel temper- 
ing" p^ , p3|]24| , |25| , p6[ , and this approach and its problems will be presented in 
Sec. 3. For the sake of contrast, Sec. 4 will then describe the 10-state Potts glass 
model. Although this model is only an abstract caricature for a real glass, it has 
the merit that quite a few results are known analytically and that Monte Carlo 
simulations are possible at T c and even at lower temperatures, if one considers 
only systems of a few hundred Potts spins. The disappointing aspect is, Sec. 4, 
that even in this very idealized case one learns relatively little about the glass 
transition of the infinite system, since one has to fight against dramatic finite 
size effects I^HI!®! 

Then we shall describe briefly (Sec. 5) a coarse-grained model of short poly- 
mer chains |p^ , pM]j3^ , |33] ,^4 35 3^,0. This beadspring model is quite successful 
in reproducing a number of experimental results qualitatively, as already exem- 
plified in Fig. [|. The cooling rates that one can reach are about 3 orders of 
magnitude smaller than for Si02 . Thus the model is very useful for testing mode 
coupling theory |3(||37]]. However, also for this system there is actually only a 
dim hope that one can get distinctly below the critical temperature T c l With 
respect to that it may be better to work with the so-called bond fluctuation 
model on the lattice " a system which allows 

to equilibrate melts at low temperatures with artificial moves @,^5 46 47|. By 
simulations of this model also the configurational entropy and its temperature 
dependence can be extracted Ej45| and thus it can be shown that the entropy 
theory of Gibbs and Di Marzio ]17fl is rather inaccurate and misleading (Sec. 6). 
Finally, an interesting variant (Sec. 7) of the bond fluctuation model will be con- 
sidered. Here one uses effective potentials that are constructed such that a real 
material is mimicked, e.g. polyethylene G§]. This trick allows that part of the 
problem of bridging the time scales is taken care of by a "time rescaling factor" 
p8| , a special translation factor between the physical time and the Monte Carlo 
time. Sec. 8 then will summarize some of the conclusions emerging from all this 
work. 



2 Towards the simulation of real glassy materials: The 
case of Si0 2 

Molten Si02 is a prototype of a network glassformer. Furthermore it is a system 
that is well suited for molecular dynamics simulations since a very well-tested 
pair potential based on quantum-chemical calculation has been developed fisf] . 
By a suitable combination of long-range Coulomb interactions and short range 
forces, chosen in the form 



V<j (r) = + A i:j exp{-B ijr ) - C y /r 6 with i, j e {Si,0}, (4) 

the effective interaction between the ions can be described reliably. Here e is the 
charge of an electron, qo — —1-2, qsi = 2.4, and the values of the parameters 
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Aij, Bij and CV,- can be found in Ref. p9| ]. This potential is able to describe 
the formation of covalent bonds without the explicit assumption of three-body 
forces, whose calculation would be very time consuming. Due to the long range 
of the electrostatic interactions, Ewald summation techniques have to be used, 
while the short range part of the potential can be cut off at a suitable radius r c . 
It turns out that r c — 5. 5 A yields good results The MD time step, however, 
must be chosen relatively small, namely St = 1.6fs. Note that the presence of 
the long range Coulombic interactions make the calculation of the forces still 
a quite CPU-intensive task. Furthermore one has to average the results over 
several independent runs in order to improve the statistics. 

In a first set of simulations, the method used to cool the sample was very 
similar to the procedure used in glass factories where the temperature is re- 
duced linearly with time t, starting from some initial temperature such that 
T(t) — Ti — jt, where 7 is the cooling rate J50| . The main difference between 
the simulation and the cooling of the real material are the actual numbers used 
here: The initial temperature that had to be chosen in the simulation was very 
high, Ti — 7000K, and cooling rates were between 7 = 10 15 K/s and 7 = 10 12 K/s 
[ pTj| . In contrast the glass factory uses typical initial temperatures of 1600K and 
cooling rates of lK/s or even less, so the simulation is at least 12 orders of mag- 
nitude off! Despite these extremely high cooling rates - which are inevitable due 
to the heavy computational burden - the generated structures are qualitatively 
reasonable. In particular one obtains random tetrahedral networks in which al- 
most all Si atom sit in the center of a tetrahedron and most of the O atoms sit 
at the corners. 

Earlier investigators (for a review see @J|,[5(|) were so bold to claim that 
such glass structures are identical to those occurring in nature, denying any 
significant dependence on cooling rate. However, as we have shown |30|, such 
a claim is foolish since one sees a pronounced dependence on cooling rate in 
many quantities, including the structure. As a typical example we show in Fig. || 
how the distribution of the length n of rings depends on 7Q. (See the figure 
caption for a definition of this length.) It is seen that over the range of 7 that is 
accessible there is a significant increase of P(n — 6) and a significant decrease 
of P(n = 3) and P(n = 4), while P(n = 5), P(n = 7) and P(n = 8) almost 
stay constant. Clearly an extrapolation of such data to the physically relevant 
cooling rate 7 = lK/s is very difficult, and perhaps not yet possible: Perhaps 
P(n = 3) and P(n = 4) are already practically equal to zero for 7 = lK/s - we 
don't really know. Even simple quantities, such as the density of the glass at low 
temperatures, are hard to predict reliably. (This problem is also complicated by 
the fact that molten Si02 has at relatively high temperatures a density anomaly 
where the thermal expansion coefficient changes sign.) 

A particular dramatic failure with extrapolations to lower values of 7 was 
encountered in an attempt to determine the cooling rate dependence of the glass 
transition temperature Tg(j). As is done in many experimental studies, one can 
fit a smooth function of temperature to the liquid branch of the enthalpy (where 
the melt has not yet fallen out of equilibrium) and another smooth function 
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to the glass branch of the enthalpy, and estimate T g [p/) from the temperature 
where these two branches intersect. Fig. ^ shows a plot of 2^(7) versus 7 - note 
the logarithmic scale for 7! - for the simulation of SiC>2. One sees, first of all, 
that there is a very strong dependence of ^(7) on 7, with T g {^) w 4000K for 
7 = 10 15 K/s, while T g {-y) has decreased down to about T g (<y) sw 2900K for 
7 = 10 13 K/s. In this range of cooling rates, the dependence of T g (^) is not linear 
in log(7_). Nonlinear variations of T s (7) that are qualitatively similar to those 
of Fig. Q have been reported in the experimental literature, too |5l|] , and are 
typically described by 



T a ( 1 )=T VF -B/[\og{ 1 A)] (5) 

This dependence can be justified by assuming that the fluid falls out of its 
(metastable) equilibrium when the time constant of the cooling, 7 , equals the 
structural relaxation time r(T) at T = T g (~f), and by using the Vogel-Fulchcr 
law from Eq. (0) for r(T), r(T) = Aexp[B/(T -T VF )}. Obviously, Eq. (|) does 
provide a very good fit to the data of the Si02 simulation, but the resulting 
Tvf = 2525K is rather unreasonable: Remember that the experimental glass 
transition temperature is T g « 1450K, the melting temperature of crystalline 
Si02 is around 2000K, and Tvf should be significantly lower than T m and even 
somewhat below T g , cf. Fig. ||. We emphasize here that the failure of the sim- 
ulation to predict T g {^f = lK/s) is not primarily due to the inaccuracy of the 
pair potential since, as will be explained in detail below, a different analysis of 
Si02 simulation data yields much more reasonable results. The failure implied 
by Fig. H simply comes from the fact that the 10-1000 picosecond timescale that 
is basically probed here is too many orders of magnitude off from the time scale 
relevant for the glass transition and that therefore an extrapolation of the results 
becomes a insecure undertaking. 

A better way to study amorphous silica is to fix density at a reasonable value, 
for instance the experimental value, and equilibrate the system at a temperature 
which is as low as possible. Present day simulations can propagate a system 
of around 8000 ions over a time span of around 20ns which allows for a full 
equilibration at T — 2750K |p4| . Longer time are accessible for smaller systems. 
However, it was found that if one has fewer than O(10 3 ) ions the results are 
plagued with finite size effects |33|. Simulating a large system over this time 
scale are on the forefront of what is feasible today, and require the use of multi- 
processor super computers such as CRAY-T3E, making use of a parallclization 
of the force calculation 1 52 5^,Q . 



This well-equilibrated melt can then be used as a starting condition for a 
cooling run at constant density. The advantage of this procedure is that a state 
at T = 2750-ftT at the correct density is much closer in local structure to the 
real glass, than the structures generated by the procedures described above, 
and hence the spurious effects of the by far too rapid quench are much less 
pronounced. This conclusion is corroborated by a comparison of the simulated 
structure factor with experiment | j55| , see Fig. |[ Given the fact that the compar- 
ison in Fig. |^ does not involve any adjustable parameter whatsoever, the agree- 
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mcnt between simulation and experiment is quite remarkable, and this reiterates 
our above conclusion that the potential used {Eq. (Q)} is accurate enough, and 
should not be blamed for discrepancies as discussed in connection with Fig. |]. 

For the temperatures at which one can equilibrate the system, i.e. here 2750K 
and higher, it is also possible to determine the self-diffusion constants of Si 
and O atoms from the simulation. This is done by calculating the mean square 
displacements (\ri(t) — T"i(0)| 2 ) = Ar 2 a (t) of the particles of type a € {Si,0}, and 
apply the Einstein relation A a (t) = 6D a t in the regime of late times where the 
dependence of A a (t) on t is in fact linear |52| , |53[|54|] } . The result is shown in Fig. [| 
where also the respective experimental data |56| , |57t are included. As one can see 
from Fig. ^, one needs to cover 16 decades, from 10~ 4 cm 2 /s to 10~ 20 cm 2 /s, 
to cover the full range including simulation results and experiments, but the 
simulation results alone are actually restricted to the first four decades of this 
range only. The straight lines fitted on this Arrhenius plot to the experiment as 
well as to the simulation show that in this case a bold extrapolation actually 
is rather successful - but of course there is no guarantee that this will work 
similarly well in other cases. 

A very interesting aspect of the temperature dependence of the diffusion 
constants is that there are strong deviations from Arrhenius behavior at very 
high temperatures. It turns out that this region is rather well described by 
a power law, D cx (T — T c ) 7 , as it is implied by mode coupling theory, see 
Eq. (||) with D cx t -1 . In fact, this conclusion is strongly corroborated by a 
detailed analysis of the intermediate scattering function <j) q (t) for wave- vector q 
and various other quantities f32| . This finding is somewhat surprising, however, 
since T c « 3330^ f5^| , i.e. far above the melting temperature of crystalline SiC^! 
Thus it is no surprise that experimental results had not given hint that mode 
coupling theory also describes a "strong" glassformer such as SiC>2 (where r 
and rj(T) follow a simple Arrhenius behavior over a wide range of temperature). 
Nevertheless, this discovery that a critical temperature exists also for SiC>2 is of 
great interest, because it suggests that the differences of the relaxation dynamics 
between different glassforming fluids are of a quantitative nature only, while 
qualitatively the behavior is always the same. 



3 Parallel tempering 

One of the major reasons for the slowing down of the dynamics of (atomistic) 
glass forming systems is that at low temperatures each atom is trapped in a 
cage formed by its surrounding neighbors. On the other hand the atom itself is 
part of a cage that trap the neighboring atoms. With decreasing temperature 
each of these cages becomes stiffer and stiffer and finally each atom can perform 
only a rattling motion, i.e. the system has become a fluid that doesn't flow 
anymore, i.e. a glass. The basic idea of the parallel tempering method is to help 
the particles to escape their local cage by supplying them with sufficient kinetic 
energy to overcome the local barrier. Originally proposed for spin models [2^p3]| , 
the method has been found to be also useful for off-lattice systems. A recent 
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review on the method can be found in Refs. p4pq |. In the following we discuss 
briefly how the method is implemented in practice. 

If we denote the Hamiltonian of the system as H = K{p) + E(q), where K 
and E are the kinetic and potential energy, respectively, and p = (pi,P2> • • • ,Vn) 
and q = (qi, q 2 , . . . , gjv) are the momenta and coordinates of the particles, we 
construct a new Hamiltonian TL as follows: 

Make M independent copies of the Hamiltonian H: Hi — K(pi) + -E(qi). 
Here the pi and qi are the momenta and coordinates belonging to the i— th 
subsystem. TL is then defined as 



M M 

H{pi, ■ ■ ■ ,p M ,qi, • - .qAf) = y^-H»(P»>q») = ^2 K (Pi) + AiE(qj). (6) 

i=l i=l 

The 1 — Ai > A 2 > . . . Am are constants which we will use later. We now make a 
molecular dynamics simulation of the Hamiltonian TL at a constant temperature 
T = fi^ 1 ■ After a certain time interval AtpT we attempt to exchange the two 
configurations m and n belonging to two neighboring systems (i.e. m — n± 1). 
Whether or not the swap of these two configurations is accepted depends on a 
Metropolis criterion with a acceptance probability 

_ / 1, Am,n < , 7 s 

Wm > n -\exp(-A m , n ), An,«>0, [n 

where A m>n — Po(A n — A m )(E(q m ) — E(q n )). Since the normal molecular dy- 
namics simulation as well as the Monte Carlo procedure on time scale At-pT fulfill 
the condition of detailed balance, the whole algorithm does so also, i.e. after a 
sufficiently long time the system composed by the subsystem will converge to a 
Boltzmann distribution. Note that in the systems with a small value of A the 
interaction between the particles is weakened (see Eq. (^)). Therefore it can be 
expected that the particles in these systems move faster than those in systems 
with a large value of A. Another way to see this is to say that each system is 
simulated at a different temperature and that periodically the temperature of 
the system is increased or decreased (hence the name of the algorithm). This 
walk in temperature space should thus allow the system to overcome the local 
barriers formed by the above mentioned cages and thus to propagate faster in 
configuration space. 

Note that this algorithm has a substantial number of parameters, all of which 
influence its efficiency considerably. In order that the acceptance probabilities 
of Eq. (0) are reasonably high, the coupling constants Ai should not be too 
different. On the other hand one wants that Am is as small as possible since this 
will lead to a fast propagation of the system at this temperature. Therefore one is 
forced to choose a relatively large value of M. This in turn is, however, not good 
for the overall performance of the algorithm since in order to be ergodic each 
configuration has to make a random walk in A— space, and the time to do this 
increases like M 2 . Last not least there is the exchange time At-pT which should 
not be too small since then the system just swaps back and forth configurations 
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that are very similar. On the other hand AtpT should also not be too large, 
since one needs these type of moves in order to explore the A— space quickly. The 
optimal choice of these parameters is currently not known and still the focus of 
research |58|. A further problem is to find out after which time the system 7i 
has really equilibrated. It seems that to guarantee this it is not sufficient that 
every subsystem has visited every point in A— space | 

should look like the one shown in Fig. ^. Furthermore we point out that it might 
be possible that a suboptimal choice of these parameters might make the whole 
algorithm rather inefficient J59| . 

If the above mentioned parameters of the algorithm are chosen well, the 
parallel tempering method can indeed speed up the equilibration of the system 
considerably. This is demonstrated in Fig. |§| where we show the mean squared 
displacement of the silicon atoms in SiC>2 as a function of time. From the figure 
we see that at the lowest temperatures the mean square displacement increases 
by about a factor of 100 faster than the corresponding curve obtained from the 
conventional molecular dynamics simulation. From the figure it becomes also 
clear that the parallel tempering slows down the dynamics of the system at high 
temperatures. This is due to the fact that these systems are coupled to the ones 
at the low temperatures and hence cannot propagate as fast anymore. 

Before we conclude this section we mention that the parallel tempering algo- 
rithm has been found to be also very efficient for the equilibration of the Potts 
glass discussed in the next section. Thus, although the algorithm might have 
some problems for certain systems or values of parameters, there are models 
where it seems to work very well. 

4 An abstract model for static and dynamic glass 
transitions: The 10-state mean field Potts glass 

In this section we are concerned with a model for which it is known exactly that 
there is a dynamical (ergodic to nonergodic) transition at a temperature To and 
a second, static, transition at a lower temperature To < Td, where a static glass 
order parameter q appears discontinuously: the infinite range p-state Potts glass 
with p > 4 |^[|]H,H,fllfl^. In this mod el, one has Potts "spin" variables 
(Tj which can take one out of p discrete values which we simply label from 1 to 
p, o~i S {1,2, ... ,p}, where i labels the "sites", i = 1, 2, . . . , N. An energy pJij 
is gained if two spins Uj, o~j are in the same state, 

i<j 

Every spin interacts with every other spin via an interaction J^- which is 
Gaussian distributed, i.e. 

exp{-(J 4 , - J ) 2 /[2(AJ) 2 }} . (9) 
Here the mean Jq and the width A J are normalized such that 



58 59]. A good random walk 



P(J K 



^(AJ) 



10 Kurt Binder et al. 



Jo = [JijU = Jo/(N - 1), (AJ) 2 = [Jl\ av - [JX = AJ/(N - 1) , (10) 

a choice that ensures a sensible thermodynamic limit. We fix the temperature 
scale byjmoosing A J = 1, and set the mean of the distribution "antiferromag- 
netic" , J = 3 — p, in order to avoid any tendency towards ferromagnetic order. 
(Note that for p = 2 this model would reduce to the standard Ising mean field 
spin glass (Sherrington-Kirkpatrick model) p7| , but we shall be concerned with 
p = 10 here.) This model, which due to the choice Eq. (^J) exhibits quenched 
random disorder already in the high temperature phase above theglass tran- 
sition, can be solved exactly in the thermodynamic limit [|60|j6l| , |6^ , |63| . |G4|j6^j6^1 . 
One finds (Fig. ||) that slightly above To the dynamic auto-correlation function 
of the spins exhibits a two-step decay, in that a plateau develops whose life-time 
diverges at To- It is important to note that this behavior is described exactly by 
mode coupling equations of the same type as they occur for the structural glass 
transition || ! This shows that this rather abstract model might be more similar 
to a real structural glass than one would expect at a first glance. At a lower 
temperature To, a static glass transition occurs I6^j6l| , |6^ , |63|j64| , |65| , |66| , where a 
static order parameter appears discontinuously. Interestingly the static response 
function does not diverge at To, i.e. the glass susceptibility is still finite here. 
The entropy does not have a jump at To, but shows only a kink. Thus there 
is no latent heat associated with this transition! A Kauzmann temperature Tr-, 
where the (extrapolated) entropy of the high temperature phase would vanish, 
also exists, but in this case clearly Tk < Tq and Tk does not have a physical 
meaning. 

It is of course interesting to know if computer simulations can identify the 
static and dynamic glass transition in a model for which one knows from the 
exact solutions p0| , |61p^ , |63] , |64[ 65 6(| that all these glass transitions do indeed 
exist. Surprisingly, the answer to this question is "no" since very strong finite 
size effects are present. In particular it is even hard to see that the above men- 
tioned plateau in the autocorrelation function develops as one approaches the 
temperature Tjj of the dynamical transition (Fig. [II]) ^8|. This is demonstrated 
in Fig. [l(] where we show the autocorrelation function of the Potts spins as a 
function of Monte Carlo time, for 160 < N < 1280. Note that this range is of 
the same order of magnitude as the particle numbers used for simulations of the 
structural glass transition, using models such as the binary Lennard-Jones fluid 
| fL2| or similar models. No evidence for strong finite size effects was ever found 
for the latter models if N was larger than s» 1000 [|68[ . Thus, a priori it is not 
at all obvious that system sizes of the order 10 3 are completely insufficient to 
characterize the dynamics of a system in the thermodynamic limit. However, 
from Fig. [l(] we must conclude that for the present system this is indeed the 
case, at least for temperatures close to the dynamical transition temperature 
To. This is in contrast with the behavior at a high temperature, e.g. T = 1.8. 
From the figure we recognize that at this temperature there are hardly any finite 
size effects and that the data have nicely converged to the thermodynamic limit 
even for modest system sizes. 
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Brangian et al. |27],^8|,^9],[30| defined a relaxation time r from the time t 
that it takes the autocorrelation function to decay to the value C(t = r) = 0.2 
(broken straight line in Fig. |o|). This time is plotted logarithmically versus 
1/T in Fig. |l|, so an Arrhenius behavior would be a straight line on this plot. 
One can see rather clearly a crossover from a power law divergence (that would 
emerge fully in the limit TV — > oo for T > Td) to the Arrhenius law at low 
T. This behavior qualitatively resembles the behavior expected for structural 
glasses where the different valleys in the rugged energy landscape for T < Tq are 
separated by finite (free) energy barriers. In contrast to this one knows that in 
the Potts glass in the limit N — > oo these barriers are truly infinite if T < To , 
and hence the dynamics is strictly nonergodic. 

Similar finite size effects affect also the behavior of static properties p7|j23 , p9 
One might wonder whether it is possible to use these finite size effects to apply 
standard finite size scaling analyses to extract reliable information on the lo- 
cation of the static transition temperature from the simulations. Unfortunately 
the answer is "no": As figure |l2| shows, the standard method J6^] of locating 
a static transition from the intersection point of the order parameter cumulant 
gives rather misleading results here since the curve seem(!) to intersect at a 
wrong temperature. Thus one must conclude that there is need to better under- 
stand finite size effects for such unconventional glass transitions as sketched in 
Fig. 0, before one can study them reliably with simulations. 
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5 The bead-spring model: A coarse-grained model for the 
study of the glass transition of polymer melts 

We now draw attention to a model which is intermediate between the abstract 
model as considered in the previous section and the chemically realistic model 
of silica melts discussed in Sec. 2. This intermediate model is a coarse-grained 
model of glassforming polymer melts. Short polymer chains are described by a 
bead-spring model, with a chain length of A^ = 10. The (effective) monomers 
interact with each other via a truncated and shifted Lennard- Jones potential, 

U LJ (r) = 4e[(a/r) 12 - (a/r) 6 } +C, r < r c = 2.2 1 ' 6 a (11) 

while Ul.j{t) = if r > r c . The constant C is chosen such that Ulj is continuous 
at r = r c . 

The spring potential present between two neighboring beads is given by 

C/fene(0 - -{k/2)Rl log[l - [l/Rof] (12) 
with the following values of the constants |u| : 

c = 1 ,a = l,k = 30,i? = 1-5. (13) 

This choice for the parameters creates frustration in the model: the minimum 
of the bond potential along the chain occurs at a position ^ m j n ps 0.97 that 
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is incompatible with the minimum position r m j n w 1.13 of the Lennard- Jones 
potential, as far as the formation of simple crystal structures is concerned. This 
conflict between these two length scales prevents crystallization very efficiently, 
and the resulting structure of the melt and the corresponding glass resembles 
corresponding experimental data very nicely, as has already been demonstrated 
in Fig. ^ 

If one carries out "slow" cooling experiments one finds that the volume per 
monomer shows at a temperature T g « 0.41 a kink |32j. This signals that the 
system has changed from the liquid branch to the glass branch and hence has 
fallen out of equilibrium. Qualitatively, the data looks again very similar to that 
of corresponding experiments ]7C[ |. However, if one compares experiment and 
simulation more quantitatively, one notes again a big disparity in the cooling 
rates: In the simulation the temperature was reduced by AT — 0.02 every 500000 
MD time steps, each time step being St = 0.002tmd with tmd = ^(m/e) 1 ^ 2 , 
m being the effective mass of the monomeric units. If one estimates that tmd 
corresponds roughly to 10 _11 s, and that T = 1 corresponds to 500K, one arrives 
at a cooling rate of AT j At « 10 9 K/s. While this estimate is three orders of 
magnitude smaller than the corresponding cooling rate for the silica melts | p0| ] , 
it is still many orders of magnitude larger than the corresponding experimental 
cooling rates. Hence also in this case there is a huge gap between the cooling 
rates accessible in simulations and those used in real experiments. 

This model yields also qualitatively very reasonable results for the relaxation 
dynamics: The self-diffusion constant can be fitted well by the Vogel-Fulcher 
law given by Eq. ([I]), with Typ ps 0.34, below the kink temperature T g ps 0.41. 
The mode coupling critical temperature is located at T c sa 0.45, above the kink 
temperature, and the ratios T c /T g and T c /Tvf are quite reasonable. Although 
in the simulation only 1200 monomers were used, a nice plateau is found in the 
intermediate incoherent scattering function (j>g(t), see Fig. [l^. Hence one can 
conclude that no strong finite size effects are present for this model. 

Also the Rouse modes JnJ which describe the mesoscopic Brownian motion 
of the polymer chains on length scales that are between monomer-monomer dis- 
tances and the coil size, are found to relax over almost two decades in T—T c with 
relaxation times that show the mode coupling power law J33J , see Fig. |l|. Only 
very close to T c , for T < 0.46, can one see small indications that the singularity 
at T c is in fact rounded off. This model has allowed many very impressive tests 
|}5|37j of mode coupling theory, similar to an often studied binary Lennard-Jones 
mixture |l2|,^2). But similar to the case in the latter model, it has so far turned 
out impossible to study temperatures for T < T c in thermal equilibrium. And 
none of these models - neither the model for Si02, nor the binary Lennard-Jones 
model [ jlH nor the present beadspring model - could provide any clarification 
about the validity of the entropy theory |l7j . 
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6 The bond fluctuation model approach to glassforming 
polymer melts 

The bond fluctuation model is an even more 

abstract model of polymers than the bead-spring model discussed in the previous 
section, since it forces the chains to "live" on a simple cubic lattice, and all 
motions on scales smaller than a lattice constant are completely suppressed. 
In this model a polymer is represented again as a chain of effective monomers 
connected by effective bonds, but now each effective monomer is described by an 
elementary cube on the lattice that blocks all 8 sites at the corners of the cube 
from further occupation (Fig. [D]). The length of the effective bonds is allowed 
to vary from 2 to \/l0 lattice constants (taken as length unit in this section). 
The only nonbonded interaction is the one of excluded volume. The dynamics 
of the random conformational changes of the real polymer is represented in a 
crude way by attempted hops of randomly chosen monomers in randomly chosen 
lattice directions. If about one half of all lattice sites are occupied, the system 
behaves like a dense melt, and even short chains with chain length N = 10 show 
already typical polymer-properties, e.g. the scaling of the radius of gyration with 
VN, etc. 

Since real polymers show with decreasing temperature an increase of the 
persistence length and hence of the chain radius, it is natural to model this effect 
by an effective potential U (I) for the length of the effective bonds, energetically 
favoring long bonds. If one chooses as a minimum of this potential U(l m i n = 3) = 
while U(l) = e — 1 for all other bond lengths I, one also incorporates "geometric 



frustration" (Fig. 15) into the model: Each bond that reaches its ground state 
wastes the four lattice sites in between the adjoining effective monomers, which 
are completely blocked for further occupation. From the point of view of packing 
as many effective monomers as possible in a dense melt on the lattice, the bonds 
that waste lattice sites are very unfavorable. Hence configurational entropy favors 
short bonds that do not waste any other lattice sites for further occupation. 
Thus a conflict between entropy and energy is created, which is responsible for 
the glass transition observed in the Monte Carlo simulations of this model. 

This model has the big technical advantage that it can be equilibrated even 
at relatively low temperatures by the so-called "slithering snake algorithm". 
In this type of Monte Carlo moves one randomly attempts to remove a bond 
from one chain end and attach it to the other chain end in a randomly chosen 
orientation [Q . Although this algorithm does not correspond to any physically 
realistic dynamics of polymers it is a perfectly admissible Monte Carlo move for 
studying equilibrium properties. Using this algorithm, thermal equilibrium can 
be established at rather low temperatures, such as T = 0.16, where after 10 7 
steps with the conventional "random hopping" algorithm the autocorrelation 
of the end-to-end vector of the chains still has not decayed below 90% of its 
starting value |4^| . If we wish to study dynamical properties of this model, we 
first perform a run with this slithering snake algorithm, to obtain initial states 
that are characteristic for thermal equilibrium. Subsequently we can start a run 
with the normal random hopping moves of the effective monomers, which thus 
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yields a physically reasonable description of the dynamics J46[ . If one estimates 
that one Monte Carlo Step per monomer corresponds to about 10~ 12 seconds in 
real time, a run of 10 7 steps would reach a physical time of 10~ 5 seconds, which 
is several orders of magnitude longer than the typical time scales accessible 
with molecular dynamics. Using this algorithm it was hence possible to make a 
very nice test of mode coupling theory [ ^2|j4^ ], resulting in T c w 0.15 while 
Tvf = 0.125 ± 0.005. However, the investigation of the relaxation dynamics in 
the regime Tvf < T < T c seems to be very difficult also in the framework of 
this lattice model, and in fact has not yet been attempted. 

Using the bond fluctuation model it was also determined how the glass tran- 
sition temperature T g depends on the length of the chain and the results are 
compatible with the law 



T g {oo)-T g (N)ocl/N. (14) 

Such a dependence has also been found experimentally |73| |, and is one of 
the most notable predictions of the entropy theory of Gibbs and Di Marzio Jl7| . 
Therefore many experimentalists believe that this theory is correct. However, 
this conclusion is premature, as a study of the configurational entropy for the 
present lattice model shows (Fig. |l6|). While the entropy does indeed decrease 
rather strongly with increasing value of inverse temperature, starting out from an 
"athermal melt" (corresponding to infinite temperature) , this decrease becomes 
slower when one approaches the vicinity of T c , and the simulation data do not 
show that the entropy vanishes, although they also cannot rule it out that this 
happens at a T far below T c . However, if one works out the Gibbs-Di Marzio 
theory |l7j] explicitly for the present lattice model (all the input parameters 
of the theory (Tt]] can also be extracted from the simulation, so there are no 
adjustable parameters whatsoever in this comparison!), one sees that the theory 
underestimates the actual entropy considerably at all temperatures. In particular 
this failure is responsible for the vanishing of the entropy at Tk ~ 0.18, which 
obviously is a spurious result, since this temperature is even higher than T c , 
well in the melt regime where the polymer system is a liquid and not a glass. In 
fact, a slightly different approximation due to Milchev |75[ renders the entropy 
nonnegative at all temperatures, but deviates now a bit from the simulation data 
in the other direction. Thus, these investigations show that although Eq. ( fuj ) 
does indeed hold it does not imply anything about the validity of the Kauzmann 
"entropy catastrophe" . 



7 Can one map coarse-grained models onto atomistically 
realistic ones? 

From the above comments it is clear that in simulations of simplified coarse- 
grained models the range of times one can span is much larger than the one for 
chemically realistic models that include atomistic detail (microseconds rather 
than nanoseconds). On the other hand, the simplified models may elucidate 
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general concepts but they fail to make quantitative predictions on the proper- 
ties of particular materials. Thus the question arises whether one can somehow 
combine the advantages of both approaches. 

An idea to do this is to make the coarse-graining process in a more system- 
atic way and to construct coarse-grained models that "remember" from which 
atomistic system they come from. For a polymer chain, coarse-graining along the 
backbone of the chain may mean that if we label the covalent bond consecutively 
(1, 2, 3, 4, 5,6,.. .) the bonds 1, 2, 3 form the effective bond /, the bonds 4, 5, 6 
form the effective bond //, etc |76]| . The potentials on the atomistic scale (e.g. 
potentials controlling the lengths of covalent bonds, the angles between them, 
the torsional angles, etc.) have then to be translated into suitable effective po- 
tentials for the length / of the effective bonds and the angle between them. 
The simplest choice would be to assume potentials of the form 

U cS {l) = ^u (l-l ) 2 ,V ee {9) = ^o(cos0-cos0 o ) 2 . (15) 

In the past potentials of this type have indeed be extracted from the proba- 
bility distributions P{1) oc e^p[-U cS (l)/k B T] , P(9) oc exp[-V cff [9)/k B T] ob- 
served in the simulations of single chains (where long range interactions need to 
be truncated, however) |76[|77|| . Of course, the effective parameters uq, Iq, vq, Oq 
are somewhat temperature dependent, and in principle one should deduce them 
from simulations of atomistically described melts containing many chains, rather 
than from single-chain simulations |7j| . The practical implementation of how one 
constructs best the effective potentials that mimic one particular material is still 
an active topic of research {4q,fT7j ■ 

A further important aspect is the question to what extent the dynamics with 
such a coarse grained system reflects the dynamics of a real chain. Here one needs 
to focus on the slowest local process, which are hops of small groups of monomers 
to a new conformation, such that a barrier of the torsional potential is crossed. 
Without such moves involving barrier crossing no conformational changes can 
occur. In a typical case, e. g. for polyethylene at T — 500K, the time scale for such 
hops is about two orders of magnitude larger than the vibration times of bond 
lengths and bond angles. Only because of this separation of time scales one can 
hope that a coarse-grained model can describe the essential features of the slow 
dynamics in the polymer melt at all, if the time units are properly rescaled. As 
shown by Tries et al. j|8|, the knowledge of the torsional potentials allows, using a 
an approach that resembles transition state theory, to construct a "time rescaling 
factor" , that gives the translation of the time unit of the Monte Carlo simulations 
(attempted Monte Carlo steps per monomer) into physical time units (Fig. |l7j ). 
One sees that for polyethylene 1 Monte Carlo step corresponds to 0.1 to lOps, in 
the temperature region of interest. At high temperatures, namely for T = 509K, 
the accuracy of the coarse-grained model of C100H202 was tested by running a 
molecular dynamics simulation of a united atom model for about a nanosecond 
(which is of the order of the Rouse relaxation time at this temperature) for 
comparison It is found that the agreement between both approaches is 
almost quantitative. The advantage of the Monte Carlo simulation of the coarse- 
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grained model is, however, that one can easily study a supercooled melt also 
at T = 250K, a temperature which is basically inaccessible to the molecular 
dynamics approach. 

If one compares the results of the coarse-grained model to experimental data, 
e. g. for the viscosity and its temperature dependence, the agreement is encour- 
agingly good but not perfect Q. One aspect which is clearly missing in the 
coarse-grained model is the description of attractive intermolecular forces. Thus, 
while this approach of mapping atomistic models to coarse-grained ones clearly 
has a great potential, there are still nontrivial problems that need to be solved. 

8 Concluding remarks 

In this brief review, the "state of the art" of computer simulations of glassy sys- 
tems was summarized. The main problem in this field is the problem of bridging 
time scales - a supercooled fluid close to the glass transitions exhibits a nontriv- 
ial dynamic behavior that extends from very fast processes (in the picosecond 
time scale range) to very slow processes (with relaxation times of the order of 
hours). Atomistic molecular dynamics simulations of chemically realistic models 
(as exemplified here for the case of molten Si02) can treat only a very small 
part of this broad range of time scales, and also special techniques such as the 
parallel tempering method can add only one or two decades to this range but not 
more. (Note also that there are still some unsolved technical problems with this 
method ]59||). While such atomistic simulations are nevertheless useful, in par- 
ticular since they complement the time range directly accessible to experiment, 
and give a very detailed insight into the interplay between structure and dy- 
namics in supercooled fluids, they clearly cannot answer questions on the nature 
of relaxation processes for temperatures close to (the experimental) T g , and the 
possible existence for an underlying static phase transition (from a metastable 
supercooled fluid to a metastable ideal glass) at a temperature Tk < T g , Also 
molecular dynamics studies of coarse-grained models for melts of short, unentan- 
gled polymer chains suffer from similar problems, although the effective cooling 
rates in these models are about a factor of 10 3 smaller than in the model for 
silica, and one can access relaxation times that are almost in the microsecond 
range. These models are very useful as a testbed for the mode coupling descrip- 
tion of the glass transition in fragile glassformers, however. Furthermore they 
have also allowed to gain very useful insight on the relaxation between the lo- 
cal motions responsible for the glass transition (cage effect etc.) and the more 
mesoscopic Brownian motion of the polymer chains (as described by the "Rouse 
modes", for instance). 

A slightly more abstract model of the same systems, the bond fluctuation 
model of glassforming polymer melts, corroborates these conclusions, although 
due to its discrete nature it is somewhat less suitable to describe the local struc- 
ture of packing effective monomers in a polymer melt or their motion on small 
scales (confined in a cage). However, this model has the merit that it allows to 
compute the temperature dependence of the configurational entropy S(T) and 
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thus to test the correctness of theories like the one of Gibbs and Di Marzio. 
While it is found that the entropy S(T) decreases significantly if the polymer 
melt approaches the glass transition, there is clear evidence that the theory of 
Gibbs and Di Marzio is quantitatively very unreliable since it underestimates 
S(T) significantly at all temperatures, and the "entropy catastrophe" that it 
predicts is clearly an artifact of inaccurate approximations. 

Finally, studies of an even more abstract model were discussed, the 10-state 
Potts glass with mean field infinite range interactions. This model has the ad- 
vantage that it is known exactly that it has a dynamical (ergodic to nonergodic) 
transition at To as well as a static transition at a (slightly) lower temperature 
To, at which a glass order parameter appears discontinuously and the entropy 
shows a kink. The conceptional disadvantage of this model, however, is that it 
has a built-in quenched random disorder (via its random exchange couplings) at 
all temperatures, unlike systems that undergo a structural glass transition, which 
have no quenched disorder in the high temperature phase (the supercooled fluid 
for T > T g ). Monte Carlo studies of this model, intended to serve as a general 
testbed for systems with both a dynamical and a static glass transition, show 
that unexpectedly large finite size effects occur, which are poorly understood. 
Thus even for this "simple" model much more work is necessary. 

While the anticipated progress in computer hardware and algorithmic im- 
provements will allow to extend the time ranges accessible in all these simula- 
tions somewhat, there is not real hope that one can bridge the desired 15 (or 
more) decades in time in this way. More promising in principle is the approach 
of providing an explicit mapping between atomistic models (which cover the 
fast processes) and coarse-grained models (which describe the somewhat slower 
processes, in the lOps to 1/is range), so that one effectively considers the same 
model system but with different approaches on different time scales. Of course, 
this idea is difficult to work out consistently in practice, and only modest first 
steps towards its realization have been taken. Much more work in this direction 
is certainly very desirable in the future. 
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Fig. 1. a) Static collective structure factor of polybutadiene at temperatures T = 4K, 
T = 16K, and T — 270K. Note that for this system the glass transition temperature 
is T g = 180A" and the critical temperature of mode coupling theory || is T c = 220K. 
The scattering background is not subtracted here, thus the zero of the ordinate axis 
is not known precisely, and the ordinate units are just measuring absolute scattering 
intensities. From Arbe et al. |l3| . b) Static collective structure factor S(q) plotted 
versus wave-vector q, for a bead-spring model of flexible polymer chains with chain 
length TV = 10. Beads interact with the potential given in Eqs. (|ri|)-(|l3|). and lengths 
are measured in units of a, temperatures in units of e. Three temperatures T = 0.2, 
0.46 and 0.52 are shown (note that T g w 0.41 and T c fa 0.45 for this model). The 
vertical lines highlight characteristic inverse length scales (related to the end-to-end 
distance R e and radius of gyration R g of the chains as well as the first maximum and 
minimum of S(q)). From Baschnagel el at. [H. 
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Fig. 2. Right figure: Schematic plot of the viscosity r/(T) of a fluid (note rj(T) oc r) vs. 
inverse temperature 1/T. The location of the melting temperature (T ra ), the critical 
temperature of mode coupling theory (T c ) HI, the glass transition temperature (T g ) and 
the Vogel-Fulcher-Kauzmann temperature [|l|Jl 5|] (To) are shown on the abscissa. The 
glass transition temperature T g is defined empirically requiring Q r/(T = T g ) = 10 13 
Poise. Two complementary concepts to explain the glass transition are indicated by the 
schematic plots on the left: The lower figure shows the time correlation function $ q (t) 
for density fluctuations at wave-vector q which according to idealized mode coupling 
theory shows at a temperature T c a nonzero "non-ergodicity parameter" f e For T 
somewhat larger than T c , <&q{i) exhibits a plateau and the "lifetime" r of this plateau 
(as well as rf) diverge as one approaches T c ||. The upper figure shows the entropy 
difference AS(T) = SWd — -^crystal, with S m = AS(T m ). The linear extrapolation 
of AS for T < T g defines the Kauzmann temperature T via AS(T = To) = fj]|. 
Adapted from Binder Q. 
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Fig. 3. Cooling rate dependence of the probability P(n) that in the network structure 
of Si02 a ring of size n is present. A ring is defined as the shortest connection of con- 
secutive Si-O elements that form a closed loop and n is the number of these segments. 
In this simulation we used 668 oxygen and 334 Si-atoms, and cooled the sample at 
constant pressure p = in an NpT simulation, cooling from the initial temperature 
T, = 7000K to the final temperature T — K. An average over 10 independent runs 
was performed, allowing to estimate the statistical errors given in the figure. From 
Vollmayr et al. JHct]. 
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Fig. 4. Effective glass transition temperature T 9 (7) plotted vs. the cooling rate 7, for 
molecular dynamics simulations of SiC>2 using the BKS potential |^] and estimating 
7^(7) from intersection points of fit functions to the enthalpy, as described in the text. 
All data are based on averages over 10 statistically independent runs. The curve is a fit 
to the function given in Eq. (^) of the text, resulting in Tvf = 2525A". From Vollmayr 
et al. lEbTl. 
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Fig. 5. Static neutron structure factor of Si02 at room temperature (T = 300K) 
plotted versus wave-vector q. The full curve is the molecular dynamics simulation of 
Ref. |52[ , using the experimental neutron scattering lengths for Si and O atoms, while 
the symbols are the neutron scattering data of Ref. [55| ]. From Horbach and Kob 
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Fig. 6. Plot of the self-diffusion constant D of silicon atoms (Si) and oxygen atoms (O) 
in molten Si02 as a function of inverse temperature. The symbols in the upper left part 
are the results from molecular dynamics simulations and the data in the lower right 
part stems from experiments j56|,[5^]. The thin straight lines show simple Arrhenius 
behavior {D tx exp(— EA/ksT)} with various choices of the activation energy Ea, 
as indicated in the figure. The vertical broken lines indicate the experimental glass 
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Fig. 7. Time dependence of the coupling constant for a parallel tempering simulation 
of liquid SiC>2. The number of particles was 336 and the number of subsystems was 32. 
Note that the shown subsystem visits all the different coupling constants several times, 
thus giving evidence that the overall system has indeed reached equilibrium. From Kob 
et al. @. 
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Fig. 8. Time dependence of the mean squared displacement of Si in Si02 at different 
temperatures. The dashed lines are from parallel tempering runs and correspond to 
temperatures 3922K, 3585K, 3235K, 3019K, and 2750K (top to bottom). The solid 
lines are from conventional molecular dynamics runs and correspond to temperatures 
6100K, 4700K, 4000K, 3580K, 3250K, and 3000K (top to bottom). From Kob et al. ||]. 
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Fig. 9. Mean-field predictions for the p-state Potts glass with p > 4. The spin glass 
order parameter, Qea, is nonzero only for T < To and jumps to zero discontinuously 
at T = To. The spin glass susceptibility xsg follows a Curie- Weiss- type relation with 
an apparent divergence at Ts < Tq. The relaxation time r diverges already at the 
dynamical transition temperature Td- This divergence is due to the occurrence of a 
long-lived plateau of height qea in the time-dependent spin autocorrelation function 
C(t). The discontinuous transition of the order parameter, however, is not accompanied 
by a latent heat. Therefore, there is no jump in the entropy at To, but only a kink occurs. 
The extrapolation of the high-temperature branch of the entropy would vanish at a 



"Kauzmann temperature" T K = [(1/4) (p - 1)/ In p] 1/2 T s w 0.988T S . From Brangian et 
al. {§. 
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Fig. 10. Time dependence of the autocorrelation function C(t) of the spins in the 10- 
state mean field Potts glass. C(t) is normalized such that C(t = 0) = 1 and C(t — > 
oo) — for T > Td- Time is measured in units of Monte Carlo steps per spin [MCS]. 
Two temperatures are shown, T = 1.8 and T = Td = 1.142 for several values of 
N. The solid horizontal line indicates the theoretical value of the Edwards- Anderson 
order parameter qeaITd) = C(t — * oo) at T — > for TV — » oo[^6[. The horizontal 
dashed line shows the value used to define the relaxation time r, C(t = r) = 0.2. From 
Brangian a/. j28[. 




Fig. 11. Arrhenius plot of the relaxation time r of the 10-state mean field Potts glass 
model for different system sizes. Error bars of r are mostly due to sample-to-sample 
fluctuations. The vertical dashed line is the location of Td where, for N — > oo, the 
relaxation time is predicted to diverge with a power-law. From Brangian et al. feSj. 
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Fig. 12. Temperature dependence of the order parameter cumulant gi(N,T) 
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1 1 + -, — %ttt — li%ll a " \ for the 10-state mean field Potts class for three choices of 
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(A r = 160, 320 and 640). (Here g is the order parameter.) The vertical straight line 

shows the location of the static transition temperature To as predicted by the exact 

solution (Si). The inset is an enlargement of the region where the three curves intersect. 

From Brangian et al. ]28[. 
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Fig. 13. Comparison of the incoherent inteAnediate scattering function (pq(t) for the 
bead-spring model at T = 0.48 and q ~ 6.9 [~ maximum of S(q), cf. Fig. ^] with various 
approximations: a Gaussian approximation (dashed line), 4>q = ex P[ — 1 2 9o(t)/6], where 
go(t) is the mean square displacement of the monomers. The mode coupling fit for 
the regime of the so-called "/3— relaxation" (solid-line) and a fit with the Kohlrausch 
function {Eq. (Q), dotted line} also are included. The non-ergodicity parameter / is 
indicated as a horizontal dashed line. From Baschnagel et al. E3. 
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Fig. 14. Variation of the relaxation time r p of the Rouse modes with the mode index 
p for the bead spring model plotted vs. T — T c , showing also a power law fit for 
p = 4(7 P = 1.83 ± 0.02). Within the error bars, this slope provides a reasonable fit for 
all p shown. From Baschnagel et al. |p3|. 



H(b) = 




Fig. 15. Sketch of a possible configuration of monomers belonging to two different 
chains in the bond fluctuation model of a polymer melt. For one monomer of the 
lower chain, an attempted move is indicated; this jump is forbidden, however, since it 
violates the excluded volume constraint. Also the choice of a two-state energy function 
is indicated, namely TL(b) = if the bondvector b equals 6 m i n = (0, 0, ± 3a) or a 
permutation thereof (a is the lattice spacing, chosen as unit of length in the following), 
and Ti(b) = e = 1 else. Note that if a bond takes a ground state bond 6 m i n it blocks 
automatically 4 sites (the 4 sites are highlighted by empty circles). From Baschnagel 
et al. §§. 
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Fig. 16. Comparison of the temperature dependence of the entropy per lattice site as 
obtained from the simulation of the bond fluctuation model (open circles) with the 
theoretical predictions of Gibbs and Di Marzio JIt]] , Flory Q| and Milchev jr^| . Note 
that the estimates for T c and Tvf are T c « 0.15 and Tvf ~ 0.125. Therefore the 
vanishing of the entropy at T » 0.18 is an artifact due to inaccurate approximations 
involved in the calculation of S(T) via the entropy theory From Wolfgardt et al. 




1/T [K" 1 ] 

Fig. 17. Temperature dependence of the time scaling factor converting the time unit 
of the Monte Carlo simulation into femtoseconds, for the case of polyethylene. The 
straight line shows that a simple Arrhenius law is a good approximation. Adapted 
from Ref. 
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